cd "/Users/tobiasbohmelt/Dropbox/Coups Spatial/PSRM Paper/2_Final Version/PSRM Replication/"

* Note: set working directory according to your location of the files.

clear
pr drop _all
set matsize 8000
set more off
version 9                                                    

program define splag_ll

args lnf mu rho sigma
tempvar A rSL
gen `A'= ones - `rho'*EIGS1
gen `rSL'=`rho'*SL1
qui replace `lnf'= ln(`A') + ln(normden($ML_y1-`rSL'-`mu', 0, `sigma'))
end

global nobs = 3497

clear
spatwmat using "inverse_coup", name(W)standardize
matrix eigenvalues eig1 imaginaryv = W
matrix eig2 = eig1'
matrix ones=J($nobs,1,1)

drop _all

use "Monadic Baseline.dta", clear

global Y effective_number
global X polity2 durable ny_gdp_pcap_kd sp_pop_totl fuel_i war_total army_centrality lag_DV country_fe2-country_fe145 year_fe2- year_fe30
global Z effective_number
mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)
matrix SL = W*Z
svmat SL, n(SL)
svmat eig2, n(EIGS)
svmat ones, n(ones)

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

ml model lf splag_ll (mu: $Y=$X) (rho:) (sigma:), technique(dfp) vce(oim)
ml init OLSb, skip
ml init rho:_cons=0 
ml init sigma:_cons=`OLSsigma'
ml max, difficult ltolerance(1e-8) tolerance(1e-8)

clear
pr drop _all
set matsize 8000
set more off
version 9                                                         

program define splag_ll

args lnf mu rho sigma
tempvar A rSL
gen `A'= ones - `rho'*EIGS1
gen `rSL'=`rho'*SL1
qui replace `lnf'= ln(`A') + ln(normden($ML_y1-`rSL'-`mu', 0, `sigma'))
end

global nobs = 3497

clear
spatwmat using "inverse_coup_lost", name(W)standardize
matrix eigenvalues eig1 imaginaryv = W
matrix eig2 = eig1'
matrix ones=J($nobs,1,1)

drop _all

use "Monadic Baseline.dta", clear

global Y effective_number
global X polity2 durable ny_gdp_pcap_kd sp_pop_totl fuel_i war_total army_centrality lag_DV country_fe2-country_fe145 year_fe2- year_fe30
global Z effective_number
mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)
matrix SL = W*Z
svmat SL, n(SL)
svmat eig2, n(EIGS)
svmat ones, n(ones)

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

ml model lf splag_ll (mu: $Y=$X) (rho:) (sigma:), technique(dfp) vce(oim)
ml init OLSb, skip
ml init rho:_cons=0
ml init sigma:_cons=`OLSsigma'
ml max, difficult ltolerance(1e-8) tolerance(1e-8)

clear
pr drop _all
set matsize 8000
set more off
version 9

program define splag_ll

args lnf mu rho sigma
tempvar A rSL
gen `A'= ones - `rho'*EIGS1
gen `rSL'=`rho'*SL1
qui replace `lnf'= ln(`A') + ln(normden($ML_y1-`rSL'-`mu', 0, `sigma'))
end

global nobs = 3497

clear
spatwmat using "Joint_Demo II.dta", name(W)standardize
matrix eigenvalues eig1 imaginaryv = W
matrix eig2 = eig1'
matrix ones=J($nobs,1,1)

drop _all

use "Monadic Baseline.dta", clear

global Y effective_number
global X polity2 durable ny_gdp_pcap_kd sp_pop_totl fuel_i war_total army_centrality lag_DV country_fe2-country_fe145 year_fe2- year_fe30
global Z effective_number
mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)
matrix SL = W*Z
svmat SL, n(SL)
svmat eig2, n(EIGS)
svmat ones, n(ones)

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

ml model lf splag_ll (mu: $Y=$X) (rho:) (sigma:), technique(dfp) vce(oim)
ml init OLSb, skip
* Initial Estimate at .0018088, but no SE. 
ml init rho:_cons=0.00001
ml init sigma:_cons=`OLSsigma'
ml max, difficult ltolerance(1e-8) tolerance(1e-8)

clear
pr drop _all
set more off
set matsize 8000
                                                     
program define splag_ll

args lnf mu rho1 rho2 sigma
tempvar A rSL1 rSL2 
gen double `rSL1'=`rho1'*SL1
gen double `rSL2'=`rho2'*SL2
scalar p1 = `rho1'
scalar p2 = `rho2'
matrix p1W1 = p1*W1
matrix p2W2 = p2*W2
matrix IpW = I_n - p1W1 - p2W2
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`rSL2'-`mu', 0, `sigma'))
end

spatwmat using "inverse_coup", name(W1)standardize

spatwmat using "Joint_Demo II.dta", name(W2)standardize

drop _all
use "Monadic Baseline.dta", clear
global nobs = 3497
matrix I_n = I($nobs)
global Y effective_number
global X polity2 durable ny_gdp_pcap_kd sp_pop_totl fuel_i war_total army_centrality lag_DV country_fe2-country_fe145 year_fe2- year_fe30
global Z effective_number
mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)
matrix SL1 = W1*Z
svmat SL1, n(SL1)
matrix SL2 = W2*Z
svmat SL2, n(SL2)

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

ml model lf splag_ll (mu: $Y=$X) (rho1:) (rho2:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
ml init rho2:_cons=0
ml init sigma:_cons=`OLSsigma'
ml max

clear
pr drop _all
set more off
set matsize 8000                                                      

program define splag_ll

args lnf mu rho1 rho2 sigma
tempvar A rSL1 rSL2 
gen double `rSL1'=`rho1'*SL1
gen double `rSL2'=`rho2'*SL2
scalar p1 = `rho1'
scalar p2 = `rho2'
matrix p1W1 = p1*W1
matrix p2W2 = p2*W2
matrix IpW = I_n - p1W1 - p2W2
qui gen double `A' = ln(det(IpW))/$nobs if _n == 1
scalar A = `A'
qui replace `lnf'= A + ln(normalden($ML_y1-`rSL1'-`rSL2'-`mu', 0, `sigma'))
end

spatwmat using "inverse_coup_lost", name(W1)standardize

spatwmat using "Joint_Demo II.dta", name(W2)standardize

drop _all
use "Monadic Baseline.dta", clear
global nobs = 3497
matrix I_n = I($nobs)
global Y effective_number
global X polity2 durable ny_gdp_pcap_kd sp_pop_totl fuel_i war_total army_centrality lag_DV country_fe2-country_fe145 year_fe2- year_fe30
global Z effective_number
mkmat $Y, matrix(Y)
mkmat $Z, matrix(Z)
matrix SL1 = W1*Z
svmat SL1, n(SL1)
matrix SL2 = W2*Z
svmat SL2, n(SL2)

qui regress $Y $X
matrix OLSb=e(b)
local OLSsigma=e(rmse)

ml model lf splag_ll (mu: $Y=$X) (rho1:) (rho2:) (sigma:)
ml init OLSb
ml init rho1:_cons=0
ml init rho2:_cons=0
ml init sigma:_cons=`OLSsigma'
ml max
